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Modeling  systems  that  are  not  inherently  isotropic,  e.g.,  extended  bilayers,  using  molecular  simu¬ 
lation  techniques  poses  a  potential  problem.  Since  these  methods  rely  on  a  finite  number  of  atoms 
and  molecules  to  describe  the  system,  periodic  boundary  conditions  are  implemented  to  avoid  edge 
effects  and  capture  long-range  electrostatic  interactions.  Systems  consisting  of  a  solvated  bilayer 
adsorbed  on  a  solid  surface  and  exposed  to  an  air/vacuum  interface  occur  in  many  experimental 
settings  and  present  some  unique  challenges  in  this  respect.  Here,  we  investigated  the  effects  of  im¬ 
plementing  different  electrostatic  boundary  conditions  on  the  structural  and  electrostatic  properties 
of  a  quartz/water/vacuum  interface  and  a  similar  quartz- supported  hydrated  lipid  bilayer  exposed  to 
vacuum.  Since  these  interfacial  systems  have  a  net  polarization,  implementing  the  standard  Ewald 
summation  with  the  conducting  boundary  condition  for  the  electrostatic  long-range  interactions  in¬ 
troduced  an  artificial  periodicity  in  the  out-of-plane  dimension.  In  particular,  abnormal  orientational 
polarizations  of  water  were  observed  with  the  conducting  boundary  condition.  Implementing  the 
Ewald  summation  technique  with  the  planar  vacuum  boundary  condition  and  calculating  electro¬ 
static  properties  compatible  with  the  implemented  electrostatic  boundary  condition  removed  these 
inconsistencies.  This  formulation  is  generally  applicable  to  similar  interfacial  systems  in  bulk  solu¬ 
tion.  ©  2011  American  Institute  of  Physics,  [doi:  10. 1063/1.3548836] 


L  INTRODUCTION 

As  many  chemical,  biological,  and  physiological  pro¬ 
cesses  sensitively  depend  on  long-range  electrostatic  forces,  it 
is  essential  that  the  molecular  descriptions  used  to  probe  these 
systems  correctly  account  for  these  forces.1  Periodic  bound¬ 
ary  conditions  (PBCs)  are  the  natural  choice  in  computer 
simulations  to  avoid  edge  effects  due  to  the  finite  number  of 
atoms  and  molecules  located  in  a  unit  cell.2  Treatments  of 
electrostatic  interactions  under  PBCs  can  strongly  influence 
calculated  properties  and,  hence,  the  physical  and  biological 
interpretation  of  the  modeled  systems.  A  number  of  different 
methods3,4  have  been  implemented  for  calculating  electro¬ 
static  interactions  in  systems  described  by  PBCs.  The  method 
of  choice,  the  Ewald  summation  method,2,3,5,6  combined 
with  the  particle  mesh  Ewald  technique,7  provides  an  efficient 
mean  to  calculate  electrostatic  interactions  for  molecular 
systems  subject  to  PBCs  and  has  been  broadly  implemented 
in  biomolecular  simulation  studies.8  Solid- supported  lipid 
bilayer  films  exposed  to  an  air/vacuum  interface  have  been 
studied  extensively  experimentally9  and  computationally.10, 11 
However,  the  implementation  of  the  Ewald  sum  in  such 
nonisotropic  systems  presents  some  unique  problems.  An 
improper  consideration  of  electrostatic  boundary  conditions 
for  such  systems  can  result  in  nonphysical  artifacts.12, 13  Here, 
we  present  an  implementation  of  the  electrostatic  boundary 
condition  that  includes  a  straightforward  correction  term  to 
remove  these  artifacts. 
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The  Ewald  summation  formula  contains  a  surface 
term  that  depends  on  the  dielectric  constant  of  the  medium 
surrounding  the  periodic  images.2,5,14,15  This  surface  term 
vanishes  for  the  conducting  boundary  condition  where  the 
dielectric  constant  of  the  surrounding  medium  is  infinite.2,5 
However,  for  vacuum  boundary  conditions  where  the  sur¬ 
rounding  medium  is  characterized  by  a  dielectric  constant 
of  1,  the  surface  term  is  expressed  as  a  function  of  the  total 
dipole  moment  or  polarization  of  the  system.2,5  The  surface 
term  has  a  slightly  different  functional  form  depending  on  the 
order  or  mode  of  summation.14, 15  It  plays  an  important  role  in 
modeling  molecular  systems  that  are  not  inherently  isotropic, 
e.g.,  extended  bilayers,  or  polarized  interfacial  systems.16  For 
systems  with  planar  interfaces,  the  planar  vacuum  boundary 
condition  where  periodic  images  are  surrounded  by  planar 
vacuum  boundaries  would  be  a  natural  choice.  For  this 
boundary  condition,  the  Ewald  summation  is  carried  out  in 
the  planar  manner,  i.e.,  summing  in  planar  directions  first 
and  then  progressing  in  the  out-of-the  plane  direction.  The 
resulting  surface  term  depends  only  on  the  total  dipole  mo¬ 
ment  in  the  out-of-the  plane  direction.  4  Yeh  and  Berkowitz16 
demonstrated  that  simulations  with  the  surface  term  cor¬ 
responding  to  the  planar  vacuum  boundary  condition  in 
conjunction  with  an  extra  vacuum  space  along  the  out-of-the 
plane  direction  reproduce  results  obtained  by  the  rigorous, 
but  computationally  expensive,  two-dimensional  Ewald 
summation  technique17  almost  exactly.  This  approximation 
holds  as  long  as  the  length  of  the  extra  vacuum  space  along 
the  out-of-the  plane  direction  is  longer  than  the  periodic  box 
lengths  in  planar  directions.16  Subsequently,  simulations  of 
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many  interfacial  systems  have  successfully  implemented  this 
surface  term  for  the  planar  vacuum  boundary  condition.18,19 
However,  the  vast  majority  of  molecular  dynamics  (MD) 
simulations  are  still  being  performed  with  the  conducting 
boundary  condition  regardless  whether  the  simulated  system 
has  a  net  polarization  or  not.  Thus,  it  is  important  to  highlight 
the  consequence  of  using  conducting  boundary  conditions  in 
simulations  of  systems  with  nonzero  net  polarization. 

It  is  important  when  simulations  are  performed  using  a 
specific  electrostatic  boundary  condition  that  the  calculations 
of  electrostatic  properties  are  consistent  with  the  implemented 
boundary  condition.  Sachs  et  al .20  and  Yeh  and  Hummer21 
introduced  an  accurate  formula  to  calculate  the  electrostatic 
potential  profile  from  simulations  performed  with  the  con¬ 
ducting  boundary  condition.  This  formula  includes  a  pre¬ 
viously  overlooked  integration  constant,  which  is  necessary 
to  enforce  the  three-dimensional  (3D)  PBC.  Interestingly,  as 
will  be  shown  below,  this  integration  constant  has  the  same 
functional  form  as  the  surface  term  for  the  planar  vacuum 
boundary  condition  but  with  an  opposite  sign.  As  a  result,  the 
traditional  formula  for  calculating  the  electrostatic  potential 
without  the  integration  constant  was  found  to  be  appropriate 
for  simulations  performed  with  the  planar  vacuum  boundary 
condition. 

Earlier  considerations  of  the  planar  vacuum  boundary  did 
not  specifically  address  the  systems  considered  here,  in  partic¬ 
ular,  those  with  intrinsically  isotropic  components  mixed  with 
a  strongly  polarized  local  component.  These  systems  with 
nonzero  net  polarization  are  encountered  even  in  the  absence 
of  charged  surfaces  or  external  electric  fields,  e.g.,  in  net- 
neutral  asymmetric  lipid  bilayers  or  quartz  crystal  microbal¬ 
ance  measurements  performed  on  polarized  crystals.22  In  this 
report,  we  examined  the  structural  and  electrostatic  properties 
of  interfacial  systems  with  nonzero  net  polarization  calculated 
with  different  electrostatic  boundary  conditions.  We  observed 
abnormal  orientational  polarizations  of  water  in  simulations 
of  a  quartz/water  interface  and  a  similar  quartz- supported  hy¬ 
drated  dimyristoylphosphatidylcholine  (DMPC)  bilayer  em¬ 
bedded  in  vacuum,  performed  while  using  the  standard  Ewald 
summation  with  the  conducting  boundary  condition.  Imple¬ 
menting  the  Ewald  summation  technique  with  the  planar 
vacuum  boundary  condition  removed  these  abnormal  polar¬ 
izations  of  water.  Here,  we  give  the  correct  formulas  for 
obtaining  electrostatic  properties  consistent  with  different 
boundary  conditions.  We  also  confirmed  that  the  modeled 
quartz- supported  hydrated  lipid  bilayer  reproduced  the  essen¬ 
tial  features  of  the  lipid  bilayer  in  bulk  water.  Finally,  we 
showed  that  the  planar  vacuum  boundary  condition  can  be  ap¬ 
plied  to  bulk  phase  simulations  of  a  quartz/water  interface  in 
order  to  remove  artifacts  introduced  by  the  conducting  bound¬ 
ary  condition. 

II.  METHODS 
A.  Simulations 

MD  simulations  were  performed  with  the  NAMD  MD  sim¬ 
ulation  program23  and  CHARMM  all  atom  force  field24  for 
DMPC.  The  TIP3P  water  model25  was  used  to  describe  water 


molecules.  The  force  field  for  a  quartz  (Oil)  crystal  was  taken 
from  the  one  developed  by  Lopes  et  al.26  A  primitive  unit  cell 
of  quartz  (Oil)  of  ~  15 -A  thickness  was  replicated  in  two  di¬ 
mensions  to  construct  a  quartz  (Oil)  crystal  with  a  surface 
area  of  36.5  x  37.0  A2,  as  in  the  previous  study  by  Lopes 
et  al 26  One  side  of  the  constructed  crystal  was  covered 
with  hydrophilic  silanols  (Si-OH),  and  silicon  atoms  on  the 
other  side  were  saturated  with  hydrogens  (Si-H).  Crystal 
atoms,  except  for  the  O-H  groups  of  silanols,  were  held 
fixed  to  maintain  the  (Oil)  crystal  geometry.  Simulations  of 
the  quartz/water  interface  and  the  quartz- supported  lipid  bi¬ 
layer  embedded  in  vacuum  were  performed  with  fixed  sim¬ 
ulation  box  sizes  using  3D  PBCs.  However,  the  box  length 
( Lz )  in  the  z-direction  normal  to  the  interface  was  chosen  to 
be  large  enough  for  the  system  to  form  interfaces  with  the 
vacuum.  This  allowed  the  systems  to  adjust  to  an  optimum 
density  without  constant  pressure  simulations.  The  penetra¬ 
tion  of  water  molecules  into  the  vacuum  phase  was  negli¬ 
gible  under  all  systems  and  conditions  studied  in  this  work. 
Short-range  interactions  outside  a  10- A  cutoff  were  truncated. 
Long-range  electrostatic  interactions  were  calculated  with  the 
particle  mesh  Ewald  method  with  and  without  the  surface 
term  for  the  planar  vacuum  boundary  condition,16  referred 
to  as  EW3DC  and  EW3D,  respectively.  The  EW3D  method 
is  equivalent  to  implementing  the  Ewald  summation  with  the 
conducting  boundary  condition.  The  contribution  to  the  po¬ 
tential  energy  (Uc)  from  the  EW3DC  correction  term  is  given 
by  the  following  surface  term  for  the  planar  vacuum  boundary 
condition: 


Uc  = 


1  2tt 
47T£o  V 


M2 


1 


Mi 


(i) 


where  £o,  V,  and  Mz  are  the  vacuum  permittivity,  the  volume 
of  the  unit  simulation  cell,  and  the  z-component  of  the  total 
dipole  moment  of  the  unit  simulation  cell.  The  contribution  to 
the  force  (Fz>i)  acting  on  atom  i  with  an  electrostatic  charge  of 
qt  from  this  term  is  as  follows: 


Fz,i  = 


dUc 

dZi 


(2) 


We  implemented  the  EW3DC  correction  term  in  NAMD  by 
adding  forces  due  to  the  following  effective  electric  field 
(£EW3DC)  z-direction16: 


£EW3DC  =  _lk'  (3) 

This  electric  field  can  also  be  related  to  the  z-component  of 
polarization  (Pz)  as  follows: 


£EW3DC  =  (4) 

£o  ’ 

Bonds  involving  hydrogens  were  constrained  with  the  SHAKE 
algorithm.27  A  time  step  of  2  fs  was  used  for  the  time 
integration.  The  temperature  was  maintained  at  310  K 
using  Langevin  dynamics  with  a  damping  coefficient  of 
10  ps-1.28  Coordinates  were  saved  at  every  1  ps  for  further 
analysis.  Errors  reported  in  this  work  were  based  on  the  95% 
confidence  interval  estimated  as  1.96  times  the  standard  error. 
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Quartz  Water 


►Z 

FIG.  1 .  A  snapshot  from  the  simulation  of  a  quartz/water  interface  embedded 
in  vacuum.  The  blue  box  represents  the  unit  simulation  cell  with  a  box  length 
in  the  z-direction  ( Lz )  of  100  A.  The  simulation  cell  is  repeated  periodically 
in  all  three  directions.  Water  is  in  contact  with  the  hydrophilic  side  of  the 
quartz  (011).  Empty  space  (vacuum)  is  used  to  separate  water  and  the  other 
side  of  the  quartz  crystal  in  the  z-direction. 


Standard  errors  were  estimated  from  an  analysis  of  five  equal 
time  blocks  from  the  trajectory  of  each  production  run. 


B.  Quartz/water  interface 

A  quartz/water  interface  was  prepared  by  adding  985  wa¬ 
ter  molecules  to  the  hydrophilic  side  of  the  crystal.  Four  sim¬ 
ulations  with  different  Lz  values  of  100,  150,  200,  and  300  A 
were  performed  with  the  EW3D  method.  The  smallest  box 
length  Lz  of  100  A  was  used  for  a  simulation  with  the 
EW3DC  method.  A  snapshot  of  the  prepared  quartz/water 
system  for  a  Lz  of  100  A  is  shown  in  Fig.  1.  A  water/vacuum 
interface  was  formed  due  to  the  extra  empty  space  between 
water  and  the  other  side  of  the  crystal  in  the  simulation  cell. 
Figure  1  shows  that  the  length  of  the  extra  empty  space  in  the 
simulation  cell  with  the  smallest  box  length  is  clearly  larger 
than  lateral  box  lengths,  which  is  required  for  the  successful 
application  of  the  EW3DC  method  to  effectively  remove  the 
periodicity  in  the  z-direction.16  Each  simulation  lasted  for 
20  ns,  of  which  the  initial  1  ns  was  treated  as  an  equilibration 
period  and  discarded. 


C.  Quartz-supported  lipid  bilayer 

A  hydrated  lipid  bilayer  composed  of  44  DMPC 
molecules  (22  molecules  on  each  leaflet)  and  1400  water 
molecules  was  prepared  in  the  simulation  cell  of  a  rectangular 


Quartz  DMPC  Water 
- -Z 


FIG.  2.  A  snapshot  from  the  simulation  of  a  hydrated  dimyristoylphos- 
phatidylcholine  (DMPC)  lipid  bilayer  deposited  on  a  quartz  crystal  surface. 
The  blue  box  represents  the  unit  simulation  cell  with  Lz  of  200  A.  The  sim¬ 
ulation  cell  is  repeated  periodically  in  all  three  directions. 


prism  with  the  CHARMM-GUI  Membrane  Builder.29  The  sur¬ 
face  area  of  the  rectangular  base  of  the  simulation  cell  was 
chosen  to  be  same  as  the  area  of  the  quartz  (Oil)  surface 
so  that  the  area  per  lipid  headgroup  was  ~60  A2,  close  to 
the  experimental  value  for  DMPC.30  The  prepared  hydrated 
DMPC  bilayer  was  equilibrated  under  NPAT  (constant  nor¬ 
mal  pressure  and  constant  lateral  surface  area  of  the  bilayer 
and  constant  temperature)  condition  with  the  CHARMM  MD 
simulation  program.31  After  a  brief  run  under  the  same  NPAT 
ensemble  with  namd,  the  lateral  periodic  boundaries  were 
gradually  changed  to  those  of  the  quartz  (Oil)  surface.  The  re¬ 
sulting  hydrated  lipid  bilayer  was  brought  closer  to  the  quartz 
(Oil)  surface,  and  the  thickness  of  the  water  layer  that  inter¬ 
faces  with  the  membrane  and  vacuum  was  doubled.  Figure  2 
shows  the  snapshot  from  MD  simulation  with  this  setup.  Both 
the  EW3D  and  EW3DC  methods  were  used  to  calculate  elec¬ 
trostatic  interactions.  The  same  Lz  of  200  A  was  used  in  both 
methods.  Figure  2  shows  that  the  length  of  the  vacuum  region 
with  Lz  of  200  A  was  clearly  larger  than  lateral  box  lengths. 
Each  simulation  lasted  for  200  ns,  of  which  the  initial  70  ns 
was  discarded  due  to  equilibration  considerations.  To  investi¬ 
gate  possible  finite  system-size  effects,  we  performed  an  addi¬ 
tional  100-ns  simulation  with  the  EW3DC  method,  where  the 
system  size  was  quadrupled  to  have  176  lipids  by  duplicating 
the  simulation  system  along  each  of  two  lateral  dimensions. 
The  value  for  Lz  was  increased  to  250  A  to  ensure  that  the 
length  of  the  empty  space  in  z-direction  is  greater  than  the 
increased  lateral  box  sizes.  Similarly,  simulation  of  a  DMPC 
bilayer  in  bulk  water  under  the  NPAT  ensemble  with  the  in¬ 
creased  system  size  was  also  performed  for  100  ns  to  compare 
properties  of  the  lipid  bilayer  in  bulk  solution  and  near  a  po¬ 
larized  quartz  interface.  The  initial  30  ns  of  these  simulations 
was  discarded  as  equilibration. 


D.  Quartz/water  interface  in  bulk  solution 

We  prepared  a  quartz/water  interface  in  bulk  solution 
without  vacuum  boundaries  by  adding  1400  water  molecules 
to  both  the  hydrophilic  and  hydrophobic  sides  of  the  quartz 
(Oil)  crystal.  In  addition  to  the  O-H  groups  on  the  hy¬ 
drophilic  side,  hydrogen  atoms  on  the  hydrophobic  side  of 
the  crystal  in  contact  with  water  were  allowed  to  move.  Sim¬ 
ulations  were  performed  for  20  ns  under  the  NPAT  ensemble. 
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The  initial  1  ns  was  discarded  as  equilibration.  We  used  both 
the  conducting  boundary  condition  with  the  EW3D  method 
and  the  planar  vacuum  boundary  condition  with  the  EW3DC 
method  to  highlight  the  physical  effects  introduced  by  the 
choice  of  boundary  condition. 

E.  Electrostatic  property  profiles  across  the  interface 


E(z)  = 


I—L-/2  Pqtfytt 

— - - h 

£o 


Pz 

£o  ’ 


(ID 


<Kz)  =  f  (z  ~  z)pq(z)dz 

£0  J-LJ. 2 


(12) 


Distributions  of  the  electrostatic  potential  and  electric 
field  along  the  z-direction  were  calculated  from  the  charge 
density  distribution  pq(z)  obtained  from  the  simulations.  Ex¬ 
pressions  for  calculating  the  electrostatic  properties  were  dif¬ 
ferent  depending  on  the  electrostatic  boundary  condition  used 
in  the  simulations.  First,  let  us  consider  the  case  when  the 
conducting  boundary  condition  was  applied  as  in  the  EW3D 
method.  Starting  from  Poisson’s  equation, 

d2<t>(z )  _  Pq{z) 

dz2  ~  s0  ’  (  ) 

where  0(z)  is  the  electrostatic  potential  at  position  z.  An  elec¬ 
tric  field  [E(z)]  was  calculated  by  integrating  Poisson’s  equa¬ 
tion  once,  as  follows: 


E(z)  = 


d<Kz) 

dz 


f—L  /2  Pq&W 

z/  - +  C 1 

£o 


(6) 


where  C\  is  an  integration  constant.  Sachs  et  al.20  and  Yeh 
and  Hummer21  showed  that  this  integration  constant  is  not 
necessarily  zero  but  can  be  determined  by  enforcing  PBCs. 
Integrating  the  above  equation  once  more,  we  obtained  the 
following  expression  for  the  electrostatic  potential: 


<Kz)  =  f  (z  ~  z')pq(z')dz' 

£0  J-LJ 2 

-Cl(z+£)+C2.  (7) 


If  we  choose  z  =  —Ez! 2  as  a  reference  point  and  set  0(— Lz!2) 
=  0,  the  integration  constant  C2  becomes  0.  Imposing  PBCs 
by  setting  4>(—Lz/2)  =  0  =  0(— Lz/2  +  Lz)  =  <fi(Lz/2)  yielded 
the  following  expression  for  C\\ 


1 

SqLz 


Pq(z')dz  ■ 


(8) 


Equation  (8)  was  expanded  to  give  the  following: 


C 1  =  - 


1 

2  £0 


(9) 


The  first  term  vanishes  because  the  net  charge  of  the  system 
is  zero.  Then, 


By  substituting  C\  =  Pz/ £ 0  into  Eqs.  (6)  and  (7),  we  obtained 
the  following  expressions: 


Similar  expressions  were  recently  obtained  by  Gurtovenko 
and  Vattulainen.32  In  their  study,  they  examined  the  elec¬ 
trostatic  potential  calculated  for  asymmetric  lipid  bilayers 
with  nonzero  net  polarization  using  the  conducting  bound¬ 
ary  condition.32  For  this  case,  they  suggested  that  calcu¬ 
lations  of  the  electrostatic  potential  should  be  done  using 
formulas  appropriate  for  the  planar  vacuum  boundary  condi¬ 
tion.  However,  Eqs.  (11)  and  (12)  should  be  used  to  calculate 
electrostatic  properties  from  simulations  performed  using  the 
EW3D  method  under  the  conducting  boundary  condition ,  re¬ 
gardless  of  whether  the  system  is  polarized  or  not,  as  we  will 
show  in  Sec.  III.  However,  when  using  the  EW3DC  method 
under  the  planar  vacuum  boundary  condition ,  the  Pz  / £q  term 
in  Eq.  (11)  is  cancelled  out  at  every  time  step  by  the  effective 
electric  field  £EW3DC  of  the  same  magnitude  but  with  the  op¬ 
posite  sign.  Therefore,  the  following  expressions  for  the  elec¬ 
tric  field  and  electrostatic  potential  should  be  used  instead: 


E(z)  = 


/-W2  Pqi-Z'W 
£0 


(13) 


4>(z)  =  — -  fZ  (z-z')Pq(z)dz.  (14) 

£0  J-LJ2 

Equations  (13)  and  (14)  are  thus  appropriate  for  the  EW3DC 
method  using  the  planar  vacuum  boundary  condition. 


III.  RESULTS  AND  DISCUSSION 
A.  Quartz/water  interface 
1.  Structural  properties 

Figure  3(a)  shows  a  comparison  of  mass  density  distri¬ 
butions  of  water  along  the  z-direction  using  different  elec¬ 
trostatic  boundary  conditions  and  box  lengths.  Overall,  these 
distributions  were  insensitive  to  changes  in  the  boundary  con¬ 
ditions  or  box  lengths  used  in  the  simulations.  In  all  cases, 
the  water  density  at  ~  10  A  away  from  the  quartz/water  inter¬ 
face  at  z  =  0  A  was  close  to  the  bulk  water  density  of  0.993 
g/cm3  at  310  K,  demonstrating  that  a  well-defined  and  sta¬ 
ble  quartz/water  interface  was  formed  with  the  addition  of  the 
water/vacuum  interface.  However,  simulations  using  smaller 
box  lengths  with  the  EW3D  method  resulted  in  slightly  lower 
water  densities  (see  Table  I).  Small  differences  in  water 
density  also  existed  near  the  water/vacuum  interface.  Lopes 
et  al  26  obtained  similar  water  density  distributions  near  the 
quartz/water  interface  from  simulations  of  water  enclosed  be¬ 
tween  quartz  surfaces.  The  distributions  are  also  similar  to 
those  obtained  by  Wander  and  Clark33  using  different  force 
field  parameters  representing  the  quartz/water  interface. 
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(a) 


with  the  EW3DC  method  but  was  larger  than  zero  for  simula¬ 
tions  with  the  EW3D  method.  These  abnormal  orientational 
polarizations  of  water  with  the  EW3D  method  were  more 
pronounced  at  the  water/vacuum  interface,  where  the  water 
density  was  low.  The  results  shown  in  Fig.  3(b)  and  Table  I 
show  that  even  though  the  difference  between  water  ori¬ 
entation  distributions  decreased  with  increasing  Lz  values 
in  the  EW3D  simulations,  noticeable  differences  persisted 
even  for  large  Lz  values.  This  indicates  that  there  were 
significant  differences  in  electrostatic  forces  experienced  by 
water  molecules  in  simulations  with  the  EW3D  and  EW3DC 
protocols. 


2.  Electrostatic  properties 

To  investigate  the  unusual  orientational  ordering  of  wa¬ 
ter  molecules  in  simulations  with  the  EW3D  method,  we  cal¬ 
culated  electrostatic  properties  across  the  quartz/water  inter¬ 
face.  Figures  4(a)  and  4(b),  show  profiles  of  the  electrostatic 
potential  and  electric  field  calculated  by  applying  Eqs.  (13) 
and  (14),  compatible  with  the  planar  vacuum  boundary  con- 


FIG.  3.  Structural  properties  of  the  quartz/water  interface,  (a)  Distributions 
of  mass  density  of  water  under  different  electrostatic  boundary  conditions, 
(b)  Average  values  of  cosine  of  the  angle  between  the  water  dipole  and  the 
z-axis  ((cos  0dipoie,z))  across  the  interface. 


The  average  cosine  of  the  angle  between  the  water  dipole 
and  the  z-axis  ((cos  $dipoie,z))  is  a  good  indicator  of  orienta¬ 
tional  polarizations  of  water  and  is  sensitive  to  the  choice  of 
electrostatic  boundary  conditions  used  in  the  simulations.13,34 
Figure  3(b)  shows  (cos  0 dipole, z)  °f  water  as  a  function  of  the 
distance  from  the  quartz  surface.  Under  isotropic  conditions, 
as  in  bulk  water  without  any  external  electric  field,  (cos 
6 dipole, z)  is  expected  to  be  zero.  In  previous  simulations 
of  water  in  contact  with  solid  surfaces  using  the  EW3DC 
method,  the  water  phase  became  isotropic  roughly  10  A  away 
from  the  surface.18,35  The  value  of  (cos  0 dipole, z)  at  the  region 
of  bulklike  water  density  (10  A  <  z  <  15  A)  was  close  to  zero 


TABLE  I.  Values  of  (cos  ^dipole, z)  and  the  average  mass  density  at  the 
bulklike  region  (10  A  <  z  <15  A)  of  water  in  the  quartz/water  interface 
as  a  function  of  different  electrostatic  boundary  conditions  and  box  size 
variations  in  the  z-direction  (Lz). 


Electrostatic 

boundary 

condition 

Lz  (A) 

(COS  ^dipole, z) 

Density  (g/cm3) 

EW3D 

100 

0.0506  (0.0001) 

0.9909  (0.0009) 

EW3D 

150 

0.0301  (0.0003) 

0.9949  (0.0008) 

EW3D 

200 

0.0214  (0.0002) 

0.9945  (0.0005) 

EW3D 

300 

0.0135  (0.0002) 

0.9967  (0.0009) 

EW3DC 

100 

0.0000  (0.0003) 

0.9972  (0.0005) 

aEW3D  and  EW3DC  refer  to  Ewald  summation  methods  without  and  with  the  correc¬ 
tion  term  for  the  planar  vacuum  boundary  condition,  respectively.  Errors  in  parentheses 
were  based  on  the  95%  confidence  interval  estimated  as  1.96  times  the  standard  error. 
Standard  errors  were  estimated  from  an  analysis  of  five  equal  time  blocks  from  the 
trajectory  of  each  production  run.  (cos  0dipoie,z)  values  are  the  average  values  of  cosine 
of  the  angle  between  the  water  dipole  and  the  z-axis. 


EW3D  lz=  100  A 
EW3DLZ=  150  A 
EW3D  Lz  =  200  A 
EW3D  Lz  =  300  A 
EW3DC  Lz=  100  A 


FIG.  4.  Electrostatic  properties  of  the  quartz/water  interface,  (a)  and  (b)  Pro¬ 
files  of  electrostatic  potential  [0(z)]  and  electric  field  [E(z)],  respectively,  cal¬ 
culated  by  applying  Eqs.  (13)  and  (14)  suitable  for  the  planar  vacuum  bound¬ 
ary  condition  to  simulation  data  generated  under  the  conducting  periodic 
boundary  condition  (EW3D).  (c)  and  (d)  The  same  profiles  are  calculated  by 
applying  Eqs.  (11)  and  (12),  i.e.,  applying  the  consistent  conducting  periodic 
boundary  condition  to  both  the  simulation  (EW3D)  and  the  calculation  of 
electrostatic  properties.  All  results  from  simulations  performed  and  analyzed 
using  the  planar  vacuum  boundary  condition  with  the  periodic  boundary  con¬ 
dition  correction  (EW3DC)  are  shown  as  solid  lines  and  were  obtained  by 
using  Eqs.  (13)  and  (14). 
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dition.  In  this  case,  the  electrostatic  potential  increased  with 
increasing  z  in  the  bulklike  water  region,  resulting  in  nega¬ 
tive  electric  field  values,  which  were  not  consistent  with  the 
positive  (cos  6 dipole, z)  values  shown  in  Fig.  3(b).  Figures  4(c) 
and  4(d),  show  the  results  obtained  by  applying  formulas  ap¬ 
propriate  to  the  electrostatic  boundary  conditions  used  in  the 
simulations,  namely,  Eqs.  (11)  and  (12)  for  simulations  us¬ 
ing  the  conducting  boundary  condition  (EW3D)  and  Eqs.  (13) 
and  (14)  for  the  simulation  using  the  planar  vacuum  bound¬ 
ary  condition  (EW3DC).  The  electrostatic  potential  profile 
obtained  with  the  EW3D  method  was  nearly  flat  in  the  re¬ 
gion  occupied  by  the  bulklike  water,  resulting  in  near  zero  but 
slightly  positive  electric  field  values,  consistent  with  the  pos¬ 
itive  (cos  0 dipole, z)  values  shown  in  Fig.  3(b).  However,  these 
relatively  weak  electrostatic  fields  were  not  strong  enough  to 
alter  the  density  distributions  of  water  molecules  embedded 
on  the  hydrogen-bonded  network  of  the  liquid  water  phase. 
In  contrast,  the  average  electric  field  in  the  region  with  the 
bulklike  water  density  with  the  EW3DC  method  was  close  to 
zero,  also  consistent  with  the  (cos  6 dipole, z)  values  shown  in 
Fig.  3(b)  and  Table  I.  This  also  demonstrates  that  the  smallest 
box  size  of  100  A  was  large  enough  to  obtain  correct  electro¬ 
static  properties  when  we  use  the  proper  correction  term  for 
the  planar  vacuum  boundary  condition,  as  implemented  in  the 
EW3DC  method. 

Gurtovenko  and  Vattulainen32  claimed  that  using 
Eq.  (12),  derived  for  conducting  boundary  conditions,  for 
simulation  systems  with  nonzero  net  dipole  moments  is  not 
appropriate  even  if  the  simulations  are  performed  with  the 
conducting  boundary  condition.  Instead,  they  suggested  using 
Eq.  (14),  which  was  derived  for  the  planar  vacuum  bound¬ 
ary  condition.  However,  our  results  in  Figs.  3  and  4  clearly 
show  that  self-consistent  electrostatic  properties  can  be  ob¬ 
tained  only  if  we  use  the  formula  corresponding  to  the  elec¬ 
trostatic  boundary  condition  used  in  the  simulation. 

B.  Hydrated  DMPC  bilayer  on  the  quartz  crystal 
1.  Structural  properties 

Figure  5(a)  shows  the  thickness  of  the  DMPC  bilayer 
represented  by  dpp,  the  average  distance  between  phosphate 
atoms  in  upper  and  lower  leaflets,  as  a  function  of  time.  Both 
10-ns  block  and  cumulative  averages  of  dpp  show  that  the  sim¬ 
ulated  systems  relaxed  within  70  ns,  our  choice  of  the  equili¬ 
bration  time.  Figure  5(b)  shows  the  contributions  to  the  mass 
density  distributions  by  DMPC  and  water  across  the  hydrated 
DMPC  bilayer  deposited  on  the  quartz  crystal.  Density  dis¬ 
tributions  calculated  from  simulations  using  either  the  EW3D 
or  EW3DC  methods  were  quite  similar.  Density  distributions 
of  water  between  the  quartz  surface  and  the  DMPC  bilayer  at 
0  A  <  z  <25  A  were  influenced  by  interactions  with  both  the 
quartz  and  the  DMPC  bilayer.  On  the  other  hand,  the  thicker 
water  layer  between  the  DMPC  bilayer  and  vacuum  phase 
contained  a  significant  bulklike  region  between  z  =  65  and 
75  A. 

Figure  5(c)  shows  the  values  of  (cos  0 dipole, z)  as  a  func¬ 
tion  of  the  distance  from  the  quartz  surface  and  demonstrates 
that  when  water  molecules  approach  the  lipid  bilayer  from  the 


FIG.  5.  Structural  properties  of  the  quartz- supported  DMPC  lipid  bilayer, 
(a)  dpp,  the  average  distance  between  phosphate  atoms  in  upper  and  lower 
leaflets,  as  a  function  of  time.  Circles  and  squares  represent  10-ns  block  av¬ 
erages  for  simulations  with  the  EW3D  and  EW3DC  methods,  respectively. 
Dashed  and  solid  lines  represent  cumulative  averages  for  the  EW3D  and 
EW3DC  methods,  respectively,  (b)  Distributions  of  mass  density  contributed 
by  the  DMPC  lipid  bilayer  and  water  molecules  calculated  from  simulations 
with  the  EW3D  (dashed  lines)  and  EW3DC  (solid  lines)  methods,  (c)  Com¬ 
parison  of  (cos  %poie, z)  values  obtained  from  simulations  with  the  EW3D 
(dashed  lines)  and  EW3DC  (solid  lines)  methods. 

bulk  phase,  the  water  dipole  initially  points  toward  the  mem¬ 
brane  at  the  interfaces  (8.4  A  <  z  <  19.3  A  or  50.2  A  <  z 
<  61.9  A).  For  water  molecules  penetrating  further  into  the 
membrane  (19.3  A  <z<  50. 2  A),  their  dipole  points  away 
from  the  membrane  center.  As  observed  for  waters  in  the  pure 
quartz/water  system,  water  molecules  in  the  simulations  with 
the  EW3D  method  were  more  polarized  compared  with  those 
with  the  EW3DC  method  even  in  the  bulklike  region. 

2 .  Electrostatic  properties 

Figure  6  shows  electrostatic  properties  calculated  for  the 
hydrated  DMPC  bilayer  supported  on  the  quartz  crystal.  As 
shown  in  Fig.  6(a),  the  electrostatic  potential  dropped  no¬ 
ticeably  across  the  membrane  region  and  the  vacuum  phase 
with  the  EW3D  method,  resulting  in  a  larger  electric  field 
[as  shown  in  Fig.  6(b)]  compared  with  those  obtained  with 
the  EW3DC  method.  At  a  first  glance,  these  potential  drops 
across  the  membrane  region  and  vacuum  region  may  seem 
physically  unreasonable.  However,  this  is  a  real  consequence 
of  having  PBCs  under  the  conducting  boundary  condition.  An 
increased  potential  in  one  part  of  the  system  (a  polarized  crys- 
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z(A) 


FIG.  6.  Electrostatic  properties  of  the  quartz-supported  DMPC  lipid  bilayer, 
(a)  and  (b)  Profiles  of  4>(z)  and  E(z)  with  the  EW3D  (dashed  lines)  and 
EW3DC  (solid  lines)  methods. 


tal  in  this  case)  has  to  be  compensated  for  elsewhere  to  satisfy 
the  imposed  PBCs. 

The  average  electric  field  near  the  region  occupied  by  the 
bulklike  water  between  DMPC  bilayer  and  the  water/vacuum 
interface  between  z  =  65  and  75  A  in  simulations  with  the 
EW3D  method  had  a  slight  positive  value  of  0.95  (±0.15) 
mV/A.  With  the  polarization  Pz  in  the  z-direction  estimated 
from  the  (cos  0 dipole, z)  value  in  the  bulklike  region  of  water  in 
Fig.  5(b),  we  can  estimate  the  dielectric  constant  e  of  water 
by  the  following  relation: 


(15) 


The  estimated  value  of  e  was  103  (±17),  in  close  agreement 
with  the  value  of  94  obtained  from  a  series  of  simulations  of 
TIP3P  bulk  water  with  different  electric  fields.36  This  con¬ 
firms  that  Eq.  (11)  produces  an  electric  field  consistent  with 
the  observed  orientational  polarization  of  water  in  simulations 
performed  with  the  EW3D  method. 
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FIG.  7.  Comparison  of  DMPC  lipid  bilayers  in  solution  (dashed  lines)  and 
with  the  quartz/water  interface  (solid  lines)  with  a  larger  system  size  having 
176  lipids.  Dotted-dashed  lines  represent  the  results  with  the  original  sys¬ 
tem  size  having  44  lipids.  Distributions  with  the  quartz/water  interface  were 
shifted  with  respect  to  the  membrane  center  to  coincide  with  the  free  sol¬ 
vated  bilayer,  (a)  Distributions  of  mass  density  contributed  by  DMPC  lipids 
and  water  molecules,  (b)  (cos  %Poie,z)  values  across  the  interface,  (c)  <f>(z) 
values,  (d)  E(z)  values,  (e)  Order  parameter  of  acyl  chain  sn-1  calculated  by 
Scd  =  l^3cos22e~l  j,  where  9  is  the  angle  between  the  C-H  vector  and  the  z- 
axis.  Circles  and  squares  represent  results  with  a  larger  system  size  for  DMPC 
bilayers  in  solution  and  with  the  quartz/water  interface,  respectively.  Pluses 
represent  the  results  with  the  original  system  size. 


C.  Comparison  of  bilayers  in  solution  and  at  the 
quartz/water  interface 

We  compared  structural  and  electrostatic  properties  of 
the  DMPC  bilayer  solvated  in  bulk  water  and  near  the 
quartz/water  interface.  Simulations  with  a  larger  system  size 
having  176  lipids  were  used  for  a  more  realistic  compari¬ 
son.  The  results  from  the  EW3DC  simulation  with  the  origi¬ 
nal  system  size  shown  in  Figs.  5  and  6  were  also  compared. 
Figure  7  shows  the  results  of  this  comparison.  Distributions 
for  the  solid- supported  bilayer  were  shifted  with  respect  to 
the  average  membrane  center  to  facilitate  comparisons  with 
the  bulk  simulation.  Figure  7(a)  shows  that  the  contributions 
to  the  mass  density  distribution  from  DMPC  were  similar 
in  bulk  water  and  near  the  quartz/water  interface.  Mass  den¬ 
sity  distributions  of  water  were  also  similar  even  though  they 
fluctuated  around  the  bulk  density  near  the  quartz  surface. 
Figures  7(b)-7(e)  show  the  comparisons  of  the  profiles  of  po¬ 
larization  of  water,  electrostatic  potential,  electric  field,  and 


lipid  acyl  chain  order  parameter,37  respectively,  in  bulk  solu¬ 
tion,  in  the  presence  of  the  quartz/water  interface,  and  with 
a  smaller  system  size.  Even  though  small  differences  ex¬ 
isted,  the  properties  across  the  bilayers  were,  overall,  very 
similar.  These  results  demonstrated  that  a  quartz- supported 
lipid  bilayer  deposited  on  top  of  a  thin  water  film  does  not 
introduce  any  gross  differences  to  a  completely  solvated  bi¬ 
layer.  Thus,  the  simulations  indicate  that  the  deposited  bi¬ 
layer  could  be  an  equivalent  model  to  study  the  properties  of 
a  lipid  bilayer  in  solution.  However,  the  properties  of  lipid  bi¬ 
layers  could  be  affected  if  a  thinner  water  layer  between  the 
solid  support  and  the  lipid  bilayer  is  modeled  as  suggested  by 
experiments9  and  simulations.11  Most  importantly,  with  the 
rigorous  treatment  of  electrostatic  boundary  conditions  de¬ 
veloped  in  this  study  for  the  quartz- supported  lipid  bilayer, 
a  variety  of  lipid/cholesterol/peptide/protein  interactions  with 
membranes  can  now  be  accurately  simulated  while  preserving 
the  inherent  anisotropy  of  the  system. 
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D.  Quartz/water  interface  in  bulk  solution 

The  planar  vacuum  boundary  condition  or  the  EW3DC 
method  is  commonly  applied  to  simulations  of  interfacial 
systems  embedded  in  vacuum  such  as  the  quartz/water 
interface  or  the  solid- supported  lipid  bilayer  described  above. 
However,  the  Pz/sq  term  in  Eq.  (11)  is  applicable  to  any  sim¬ 
ulations  performed  under  conducting  boundary  conditions, 
regardless  whether  the  simulated  system  is  in  bulk  solution 
or  embedded  in  vacuum.  Therefore,  we  investigated  whether 
artificial  PBC-induced  orientational  ordering  was  present  for 
a  fully  solvated  quartz  interface  using  conducting  boundary 
conditions  and  how  it  can  be  removed  by  implementing  the 
EW3DC  planar  vacuum  boundary  condition.  Figures  8(a) 
and  8(b),  show  a  comparison  of  distributions  of  mass  den¬ 
sity  and  orientational  polarization  of  water,  respectively, 
across  the  quartz/water  interface  in  bulk  solution  obtained 
from  simulations  with  the  EW3D  and  EW3DC  methods. 
Density  distributions  of  water  near  the  hydrophilic  side  (z 
>  7  A)  of  the  quartz  crystal  for  both  methods  were  similar 
to  those  shown  in  Fig.  3(a)  for  the  quartz/water  interface 
embedded  in  vacuum.  However,  the  unusual  ordering  of 
water  with  the  EW3D  method  using  the  conducting  boundary 


FIG.  8.  Structural  and  electrostatic  properties  of  the  quartz/water  interface 
in  bulk  solution  without  vacuum  interfaces.  Results  with  the  EW3D  and 
EW3DC  methods  are  shown  by  the  dashed  and  solid  lines,  respectively.  Dis¬ 
tributions  were  calculated  with  respect  to  the  center  of  the  quartz  crystal, 
(a)  Distributions  of  mass  density  contributed  by  water  molecules,  (b)  (cos 
^dipole, z)  values  across  the  interface,  (c)  tf>(z)  values,  (d)  E(z)  values.  The  inset 
in  C  shows  a  magnified  view  of  (f>(z)  profiles  in  the  bulklike  region  of  the 
water  phase. 


condition  in  Fig.  8(b)  was  more  pronounced  compared 
with  when  embedded  in  vacuum.  In  addition,  the  water 
density  distribution  near  the  hydrophobic  side  (z  <  —  7  A) 
in  Fig.  8(a)  was  also  affected  by  the  choice  of  the  electrostatic 
boundary  condition  due  to  the  high  degree  of  ordering  of 
water  with  the  EW3D  method.  The  distribution  of  water 
density  near  the  hydrophobic  side  was  also  slightly  different 
compared  with  that  obtained  from  a  previous  simulation26 
where  electrostatic  interactions  were  truncated.  As  expected, 
the  unusual  ordering  of  water  was  completely  removed  when 
we  used  the  planar  vacuum  boundary  condition  with  the 
EW3DC  method.  Distributions  of  the  electrostatic  potential 
and  electric  field  shown  in  Figs.  8(c)  and  8(d),  calculated  with 
formulas  corresponding  to  the  employed  boundary  condition, 
were  consistent  with  the  (cos  6 dipole, z)  values  shown  in 
Fig.  8(b).  A  similar  relationship  between  (cos  0 dipole, z)  and 
the  electric  field  was  observed  in  the  previous  simulation  of 
TIP3P  bulk  water.36  Our  results  demonstrate  that  the  planar 
vacuum  boundary  condition  can  be  applied  to  more  general 
interfacial  systems  in  bulk  solution  to  remove  the  nonphysical 
ordering  observed  under  the  conducting  boundary  condition. 

IV.  SUMMARY  AND  CONCLUSIONS 

The  general  system  consisting  of  a  solvated  bilayer  ad¬ 
sorbed  on  a  solid  surface  and  exposed  to  an  air/vacuum  inter¬ 
face  occurs  in  many  experimental  settings  and  presents  some 
unique  computational  challenges.  In  this  study,  we  examined 
how  the  structural  and  electrostatic  properties  of  polarized  in¬ 
terfacial  systems  are  affected  by  the  choice  of  the  electrostatic 
boundary  condition  and  the  importance  of  analyzing  electro¬ 
static  properties  with  the  correct  formulas  corresponding  to 
the  deployed  boundary  conditions.  We  showed  that  imple¬ 
menting  the  standard  Ewald  summation  with  the  conducting 
boundary  condition  for  interfacial  systems  with  a  net  polar¬ 
ization  introduced  an  abnormal  orientational  polarization  of 
water.  Our  analyses  suggest  that  the  polarized  component  in 
the  simulated  system  (e.g.,  a  polarized  quartz  crystal  or  an 
asymmetric  bilayer)  affects  the  properties  of  other  parts  of 
the  system  by  a  compensation  mechanism  when  we  apply  the 
conducting  boundary  condition.  Implementing  the  Ewald 
summation  technique  with  the  planar  vacuum  boundary  con¬ 
dition  and  evaluating  electrostatic  properties  with  formulas 
corresponding  to  the  electrostatic  boundary  condition  does 
not  introduce  any  such  compensatory  effect,  and  the  influence 
of  the  polarized  component  is  more  realistically  accounted 
for.  We  also  confirmed  that  the  modeled  quartz- supported  hy¬ 
drated  lipid  bilayer  reproduced  the  essential  features  of  the 
lipid  bilayer  in  bulk  water. 
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